#峰值提取函数
#基于本段函数完成数据中所有峰值的提取
import numpy as np
import matplotlib.pyplot as plt
def RS_Peak_abstract(frequencys,amplitudes,ratio):
    if len(frequencys)!=len(amplitudes):
        print("当前数据错误，频点与幅值不对应，不能进行包络提取")
    plt.figure(1)
    plt.loglog(frequencys,amplitudes)
    plt.rcParams['font.sans-serif'] = ['SimHei'] 
    plt.xlabel("Frequency/MHz")
    plt.ylabel("Amplitude/mA")
    plt.title("原始数据")
    #首先进行曲线的平滑，平滑窗口为4%
    # for i in range(len(frequencys)):
    #     amplitudes_sum=[]
    #     for j in range(len(frequencys)):
    #         if frequencys[j]>frequencys[i]*0.9 and frequencys[j]<frequencys[i]*1.1:
    #             amplitudes_sum.append(amplitudes[j])
    #     amplitudes[i]=sum(amplitudes_sum)/len(amplitudes_sum)
    # plt.figure(2)
    # plt.plot(frequencys,amplitudes)
    # plt.rcParams['font.sans-serif'] = ['SimHei'] 
    # plt.xlabel("Frequency/MHz")
    # plt.ylabel("Amplitude/mA")
    # plt.title("平滑数据")
    # plt.show()
    
    
        
    #开始提取峰值
    amplitudes_peak=[]
    frequencys_peak=[]
    for i in range(1,len(frequencys)-1):
        if amplitudes[i]>amplitudes[i-1] and amplitudes[i]>amplitudes[i+1]:
            print(f"当前已经找到峰值频率{frequencys[i]},此时的峰值为{amplitudes[i]}")
            frequencys_peak.append(frequencys[i])
            amplitudes_peak.append(amplitudes[i])
    plt.figure(2)
    plt.loglog(frequencys_peak,amplitudes_peak)
    plt.rcParams['font.sans-serif'] = ['SimHei'] 
    plt.xlabel("Frequency/MHz")
    plt.ylabel("Amplitude/mA")
    plt.title("当前峰值提取结果")
  

    #基于当前峰值开始电平的设置
    #原始数据是frequencys和amplitude，峰值数据是frequencys_peak,amplitudes_peak
    for i in range(len(frequencys)):
        for j in range(len(frequencys_peak)):
            if frequencys[i]>frequencys_peak[j]*(1*ratio) and frequencys[i]<frequencys_peak[j]*(1/ratio):
                if j==0:
                    amplitudes[i]=amplitudes_peak[j]
                elif amplitudes_peak[j]>amplitudes_peak[j-1]:
                    amplitudes[i]=amplitudes_peak[j]
                else :
                    amplitudes[i]=amplitudes_peak[j-1]

    plt.figure(3)
    plt.loglog(frequencys,amplitudes)
    plt.rcParams['font.sans-serif'] = ['SimHei'] 
    plt.xlabel("Frequency/MHz")
    plt.ylabel("Amplitude/mA")
    plt.title("当前标准限值")
    plt.show()
    

    


if __name__=="__main__":
    x=np.loadtxt(r"E:\pythonproject\simulation\HIRF\传导敏感测试数据.txt")
    #数据处理，原数据是dbuA，当前处理为mA，如果未来数据是mA，直接略过此步
    frequencys=x[:,0]
    amplitudes=x[:,1]
    amplitudes=[10**((amplitudes[i]-60)/20) for i in range(len(amplitudes))]
    combined =zip(frequencys,amplitudes)
    file = open(r"E:\pythonproject\simulation\HIRF\data_llsc.txt","w")
    for data in combined:
        file.write(f"{data[0]} {data[1]}\n")
    file.close()
    RS_Peak_abstract(frequencys,amplitudes,0.5)
    
    

    
    
    

        





